%stns = [1011,1012,2004,2010,2013,2016,2018,2020,3005,3006,3010,3014,3015,4001,4002,4003,4015,4017,4019];
minlength = 1025;

% MANOJ OFFICE
% path1 = 'C:\Manoj\ngri-dell01\data\sel_ts\xy\';
% pathout = 'C:\Manoj\PhD\final\xy-rb\';
% [A,B] = xlsread('c:\manoj\Phd\finalxyts.xls');

% MANOJ HOME
path1 = 'C:\Manoj\ngri-dell01\data\ts\';
pathout = 'C:\Manoj\phd\finalxy_ohneANN\';
[A,B] = xlsread('c:\manoj\phd\finalxyts.xls');

L = A(:,5) ~= 0;
A(L,:) = [];
B(L,:) = [];
 
stns = union(A(:,1),A(:,1))'; %eleminates duplicates


for istations = 1:length(stns),
%   for isessions = 1:max(A(:,3)),
        L = A(:,1) == stns(istations);% & A(:,3) == isessions;
        if isempty(B(L,4)) ~= 1,
            files = char(B(L,4));
              for nfiles = 1:sum(L),
                  fid = fopen([path1 files(nfiles,:)],'rb');
                 if fid == -1,
                     disp('Error in readind');
                     disp([path1 files(nfiles,:)]);
                    break;
                    end;

              file = files(nfiles,:);
disp('Rts-in -> reading time series');
[TS,fid] = rts_in(fid); % reading time series file (rts format - metronix)
fclose(fid);
if TS.err == -1,
    disp('error in reading a ts file');
    break;
end;

if length(TS.matrix) < minlength,
    fprintf('Insufficient data points (%d) skipping file %s\n',length(TS.matrix),file);
    break;
end;
ProcDef = SetDefault(TS.head.band); % setting defaults. ProcDef is a structured array
ProcDef.file1=file;
ProcDef.file2=file;
ProcDef.len = floor(length(TS.matrix)/ProcDef.block);
ProcDef.SelStacks = ones([1 ProcDef.len]);  
fprintf('File %s, Band %d, Stacks %d Nfft %d Overlap %d\n', ProcDef.file1,ProcDef.nband,ProcDef.len,ProcDef.block,...
ProcDef.overlap);
fprintf('Processing parameters: Section Averaging %s, Band Averaging %s\n',ProcDef.Method,ProcDef.BandAv);

if ProcDef.nband == 4 & ProcDef.DoANN == 1,
	ProcDef = SelANN(TS,ProcDef);
end;                                               

if ProcDef.BandAv == 'R', %robust band averaging 
   [SPMatrix] = GetRobSpectra(TS,ProcDef);
else, % Parzen band averagibg
   [SPMatrix] = GetSpectra(TS,ProcDef);
end;

temp = size(SPMatrix); % This is required, as number of overlapped ts sections may not be just 2*len-1
ProcDef.len = temp(2);
ProcDef.Weights = ones([ProcDef.nfrq ProcDef.len]); % Weights are modified in robust averaging routine(s)
%SPMatrix = ScaleSPM(SPMatrix,ProcDef);
disp('NormSPM-> Normalize the SPMs wrt power in horz. mag fields');
SPMatrix = NormSPM(SPMatrix,ProcDef); % Normalize the SPMs wrt power in horz. mag fields
switch ProcDef.Method,
case 'LS', % Least Square
   disp('AverageSPM (Non robust)-> Averaging  spectra sets....');
   SPM2  = AverageSPM(SPMatrix,ProcDef);
	disp('TransFn   - > Computing impedances....');
   tf2=TransFn(SPM2,ProcDef);
case 'RB', % Robust M
   [tf2,ProcDef,SPM2] = robspm(SPMatrix,ProcDef);
case 'CT', % Coherency Thrashold
   [SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,2);
   [tf2]=tf_dn(SPM2);
   disp('Warning Errors are not computed\n');
case 'CO', % Preferential Coherency Stacking
	[SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,1);
   [tf2]=tf_dn(SPM2);
   disp('Warning Errors are not computed\n');
case 'JK', % Jack knife impedance & error 
   [tf2,dof]=jknf(SPMatrix,ProcDef);
end;

CollectINFO(nfiles).bands = ProcDef.nband;
CollectINFO(nfiles).files = ProcDef.file1;
CollectINFO(nfiles).method = ProcDef.Method;
CollectINFO(nfiles).BandAv = ProcDef.BandAv;
CollectINFO(nfiles).Nfft = ProcDef.block;
CollectINFO(nfiles).ANN = ProcDef.DoANN;
CollectINFO(nfiles).overlap = ProcDef.overlap;

    if exist('FSPM') == 1,
        FSPM=addspm(FSPM,SPM2);
        else,
        FSPM=SPM2;
       end;
        clear SPM2 SPMatrix TS tf2 dof;
end;


FSPM = sortspm(FSPM);



if ProcDef.WriteEdi == 1,
   fprintf('Edi-out -> writing out edi\n');
   ediheader = ['C:\Manoj\ngri-dell01\data\edihead\' 'VP' num2str(stns(istations)) '.txt'];
  % ediheader = ['C:\Manoj\data\edihead\' 'VP' num2str(stns(istations)) '.txt'];
   fkk = fopen(ediheader,'rt');
   if fkk == -1,
       fprintf('Unable to open edi header file: %s error in edi_out\n',ediheader);
   end;
% file1 = ['VP' num2str(stns(istations)) num2str(isessions) '.edi'];
  file1 = ['VP' num2str(stns(istations)) '.edi'];
   fid = fopen([pathout file1],'wt');
   if fid == -1,
      fprintf('Unable to open %s error in edi_out header from %s\n',[pathout file1],ediheader);
      return;
   end;
   while feof(fkk) == 0,
       S = fgets(fkk);
       if findstr(S,'>INFO') == 1,
           fgets(fkk);fgets(fkk);fgets(fkk);fgets(fkk);
           fprintf(fid,'%s',S);
           fprintf(fid,'  EDI file from MTPROCN a matlab software for MT TS processing\n');
           fprintf(fid,'  Processing options\n');
           fprintf(fid,'  Files     Band Nfft Method overlap ANN\n');
           for iii = 1:length(CollectINFO),
               fprintf(fid,'%20s %5d %5d %5s %5d %5d\n',CollectINFO(iii).files,CollectINFO(iii).bands,CollectINFO(iii).Nfft,CollectINFO(iii).method,CollectINFO(iii).overlap,CollectINFO(iii).ANN);
           end;
           fprintf(fid,'Processed on %s\n\n',date);
       else,
       fprintf(fid,'%s',S);
   end;
   end;
   fprintf(fid,'\n');
   FSPM.SiteID = ['VP' num2str(stns(istations))];
   [fid,er]=edi_out(FSPM,fid);
   fprintf('\n\nWrote EDI file %s\n\n',file1);
   fclose(fid);
   fclose(fkk);
end;
clear FSPM ProcDef CollectINFO;
end;
%end;
end;